#pragma once
#include "global.h"
#include <bits/stdc++.h>
//#include <eigen3/Eigen/Eigen>

using namespace std;
//using namespace Eigen;

//int solveEqu(void);
//int updateInit(void);

/*
int solveEquTemp(void){
    //solve
    {
        Matrix<double,3,3> A;
        Vector3d b;
        for(int i=0;i<3;i++){
            for(int j=0;j<3;j++){
                A(i,j)=equ.matA[i][j];            
            }
        }
        for(int i=0;i<3;i++){
            b(i)=equ.vecB[i];
        }
        Vector3d x;
        x=A.lu().solve(b);
        //cout<<x<<endl;
        for(int i=0;i<equ.M;i++){
            equ.lambda[i]=x(i);
        }
        equ.W=x(equ.M);
        for(int i=equ.M+1;i<equ.M+equ.L+1;i++){
            equ.X[i-equ.M-1]=x(i);
            if(equ.X[i-equ.M-1]<0) equ.X[i-equ.M-1]=eps;
        }
    }
    //calc n
    {
        equ.XG=equ.cG*equ.W;
        cout<<"X_g="<<endl<<"\t"<<equ.XG<<endl;
        for(int j=equ.L;j<equ.N;j++){
            equ.X[j]=-equ.Phi[j]+equ.c[j]/equ.cG*equ.XG;
            for(int k=0;k<equ.M;k++){
                equ.X[j]+=equ.c[j]*equ.A[k][j]*equ.lambda[k];
            }
            if(equ.X[j]<0) equ.X[j]=eps*eps;
        }
        //cout
        cout<<"X="<<endl;
        for(int i=0;i<equ.N;i++){
            cout<<"\t"<<equ.X[i];
        }
        cout<<endl;
    }
    return 0;
}
*/

int solveEqu(){
    int n;
    double A[maxCompoundNum+maxEleNum+1][maxCompoundNum+maxEleNum+1];
    double L[maxCompoundNum+maxEleNum+1][maxCompoundNum+maxEleNum+1];
    double U[maxCompoundNum+maxEleNum+1][maxCompoundNum+maxEleNum+1];
    double b[maxCompoundNum+maxEleNum+1];
    double x[maxCompoundNum+maxEleNum+1];
    double y[maxCompoundNum+maxEleNum+1];
    //init
    
    n=equ.M+equ.L+1;
    for(int i=0;i<n;i++){
        for(int j=0;j<n;j++){
            A[i][j]=equ.matA[i][j];
        }
        b[i]=equ.vecB[i];
    }
    
    /*
    n=3;
    for(int i=0;i<3;i++) A[0][i]=A[i][0]=i+1;
    A[1][1]=A[2][2]=5;
    A[1][2]=2;A[2][1]=1;
    b[0]=14;b[1]=18;b[2]=20;
    */
    //calc diag of L
    for(int i=0;i<n;i++){
        L[i][i]=1;
    }
    //calc first row of U
    for(int i=0;i<n;i++){
        U[0][i]=A[0][i];
    }
    //calc first column of L
    for(int i=1;i<n;i++){
        L[i][0]=A[i][0]/U[0][0];
    }
    //calc U and L
    for(int r=1;r<n;r++){
        for(int i=r;i<n;i++){
            U[r][i]=A[r][i];
            for(int k=0;k<=r-1;k++){
                U[r][i]-=L[r][k]*U[k][i];
            }
            if(r==n||i==r) continue;
            L[i][r]=A[i][r];
            for(int k=0;k<=r-1;k++){
                L[i][r]-=L[i][k]*U[k][r];
            }
            L[i][r]/=U[r][r];
        }
    }
    //cout
    cout<<"U="<<endl;
    for(int i=0;i<n;i++){
        for(int j=0;j<n;j++){
            cout<<"\t"<<U[i][j];
        }
        cout<<endl;
    }
    cout<<"L="<<endl;
    for(int i=0;i<n;i++){
        for(int j=0;j<n;j++){
            cout<<"\t"<<L[i][j];
        }
        cout<<endl;
    }
    //calc y
    y[0]=b[0];
    for(int i=1;i<n;i++){
        y[i]=b[i];
        for(int k=0;k<i;k++){
            y[i]-=L[i][k]*y[k];
        }
    }
    //calc x
    x[n-1]=y[n-1]/U[n-1][n-1];
    for(int i=n-2;i>-1;i--){
        x[i]=y[i];
        for(int k=i+1;k<n;k++){
            x[i]-=U[i][k]*x[k];
        }
        x[i]/=U[i][i];
    }
    //cout
    cout<<"y="<<endl;
    for(int i=0;i<n;i++)
        cout<<"\t"<<y[i];
    cout<<endl;
    cout<<"x="<<endl;
    for(int i=0;i<n;i++)
        cout<<"\t"<<x[i];
    cout<<endl;
    //trans
    //trans
    for(int i=0;i<equ.M;i++){
        equ.lambda[i]=x[i];
    }
    equ.W=x[equ.M];
    for(int i=equ.M+1;i<equ.M+equ.L+1;i++){
        equ.X[i-equ.M-1]=x[i];
        if(equ.X[i-equ.M-1]<0) equ.X[i-equ.M-1]=eps;
    }
    {
        equ.XG=equ.cG*equ.W;
        cout<<"X_g="<<endl<<"\t"<<equ.XG<<endl;
        for(int j=equ.L;j<equ.N;j++){
            equ.X[j]=-equ.Phi[j]+equ.c[j]/equ.cG*equ.XG;
            for(int k=0;k<equ.M;k++){
                equ.X[j]+=equ.c[j]*equ.A[k][j]*equ.lambda[k];
            }
            if(equ.X[j]<0) equ.X[j]=eps*eps;
        }
        //cout
        cout<<"X="<<endl;
        for(int i=0;i<equ.N;i++){
            cout<<"\t"<<equ.X[i];
        }
        cout<<endl;
    }
    return 0;
}

int updateInit(void)
{
    for(int i=0;i<equ.N;i++){
        equ.c[i]=equ.X[i];
    }
    return 0;
}
